using namespace std;
//-------------

bool v,v1; // parametr wywolania (pokazuje co robi program)
bool m; string nazwam="macierz.log"; // parametr wywolania (zapisuje macierz do pliku)
string nazwap="przemieszczenia.rlt"; // plik wynikowy dla przemieszczen
string nazwao="odksztalcenia.rlt";
string nazwan="naprezenia.rlt";
string nazwaoo="wynik";
string nazwa="pusty";
string nazwa_msh="pusty";
string nazwa_bc="pusty";
string nazwa_m="pusty";
//------------------------ STRUCTY --------------------------
struct wezel_
{
	double x,y;
        int id;
};

struct koneksja_
{
	string typ_el;//typel-typ elementu
	int nr, wez_tw1, wez_tw2 ,wez_tw3, obszar;

 
};

struct brzeg_
{
	int nr, il_wez,ob;
	int * wt; // wezly tworzace
	bool utw_x,
             utw_y,
             obc_x,
             obc_y;
	double wobc_x,
               wobc_y, 
               wutw_x,
               wutw_y;
};

struct obszar_
{	
	string opis;
	int nr,il_el;
	int * el_tw;
};

struct bc_
{
	int nr_brzegu,
            nr_obszaru,
            il_wezlow,
            il_warunkow;
};

struct odksztalcenie_
{
   double Ex,Ey,Ez, Qxy;
};

struct naprezenie_
{
   double Sx,Sy,Sz, Yxy;
};
//---------

class Siatka
{
private:
	string opis;
	int liczba_wymiarow, 
	    liczba_wezlow, 
	    liczba_elementow, 
	    liczba_brzegow, 
	    liczba_obszarow;
	wezel_ * wspolrzedne;
	koneksja_ * koneksje;
	brzeg_ * brzegi;
	obszar_ * obszary;
	double w801; // modul younga
	double w802; // modul wzmocnienia
	double w803; // gran plast
	double w804; // liczba poissona
	double w805; // 
	double f1,f2,f3;
	bool psn,pso;
        double **mac;
       	bc_ *bc;     
public:
	void wczytaj_msh(string nazwa_pliku);
	void wczytaj_bc(string nazwa_pliku);
	void wczytaj_m(string nazwa_pliku);
	

	void setPSN() { psn=1; pso=0; }
	void setPSO() {	pso=1; psn=0; }
	
	double w_801_() { double temp=w801; return temp; }
	double w_804_() { double temp=w804; return temp; }
	
	void wspol_sprez(double e,double v);
	double * wspolczynniki(int nrele);
	void Macierz();
	void Obciazenia();
	void Utwierdzenia();
	double *Gauss();
        odksztalcenie_ * Odksztalcenia(double v);
        naprezenie_ * Naprezenia (double v, odksztalcenie_ *odk);
        void ZapiszOdksztalcenia (odksztalcenie_ *o);
        void ZapiszNaprezenia (naprezenie_ *n);
        void ZapiszPrzemieszczenia (double *t);


};
//-----------------------------------------------------------

//============= Wczytywanie MSH ===============

void Siatka::wczytaj_msh(string nazwa_pliku)
{
  nazwa_pliku+=".msh";
  
  if (v==1) cout << "Wczytywanie pliku msh:\n==============================" <<endl;

  string tmps;
  int tmpi;

  ifstream plik;
  plik.open(nazwa_pliku.c_str());
  
  plik>>tmps;
  if(tmps!="<SIATKA>")
  { cerr<<"Brak znacznika <SIATKA> w msh"<<endl;  exit(1);	}
  if (v==1) cout << "<SIATKA>"<<endl;
  plik>>tmps;
  opis="";
  if(tmps!="[opis]"&&tmps!="[OPIS]")
  { cerr<<"Brak znacznika [opis]) w msh"<<endl;	exit(1); }
  else
	{
          if (v==1) cout << "[opis]";
          plik >> tmps;
          while (tmps!="[liczba" && tmps!="[LICZBA") 
             {
                opis+=" "+tmps;
                plik >> tmps;
             }
         if (v==1) cout << opis <<endl;
        }
	
  plik >> tmps;
  if(tmps!="wymiarow]" && tmps!="WYMIAROW]")
  {
   cerr<<"Brak znacznika [liczba wymiarow] w msh"<<endl; exit(1);	
  }
   else 
  {
   plik>>tmpi;	if(tmpi!=2)
     {	cerr<<"Nieprawidlowa Liczba Wymiarow w msh"<<endl; exit(1); }
   liczba_wymiarow=tmpi;
   if (v==1) cout << "[liczba wymiarow]\n"<< liczba_wymiarow<<endl;
  }

  plik>>tmps;
  if(tmps!="[liczba" && tmps!="[LICZBA")
  { cerr<<"Brak znacznika [liczba wezlow] w msh)"<<endl; exit(1); }
   
  plik >> tmps;
  if(tmps!="wezlow]" && tmps!="WEZLOW]")
  { cerr<<"Brak znacznika [liczba wezlow] w msh)"<<endl; exit(1); }
  else	{   plik>>tmpi;	liczba_wezlow=tmpi; 
  if (v==1) cout << "[liczba wezlow]\n"<< liczba_wezlow <<endl;	}

  plik>>tmps;
  if(tmps!="[liczba" && tmps!="[LICZBA")
  { cerr<<"Brak znacznika [liczba elementow] w msh"<<endl; exit(1);}
  plik>>tmps;
  if(tmps!="elementow]" && tmps!="ELEMENTOW]")
  {cerr<<"Brak znacznika [liczba elementow] w msh"<<endl;  exit(1);}
  else
   {
     plik>>tmpi;     liczba_elementow=tmpi; 
 if (v==1) cout << "[liczba elementow]\n"<< liczba_elementow<<endl;
	}

  plik>>tmps;
  if(tmps!="[liczba" && tmps!="[LICZBA")
  { cerr<<"Brak znacznika [liczba brzegow] w msh "<<endl; exit(1); }
   
  plik>>tmps;
  if(tmps!="brzegow]" && tmps!="BRZEGOW]")
  {  cerr<<"Brak znacznika [liczba brzegow] w msh"<<endl; exit(1); }
  else
   {  plik>>tmpi;	liczba_brzegow=tmpi;	
 if (v==1) cout << "[liczba brzegow]\n"<< liczba_brzegow<<endl;
   }

  plik>>tmps;
  if(tmps!="[liczba" && tmps!="[LICZBA")
  { cerr<<"Brak znacznika [liczba obszarow] w msh"<<endl;  exit(1);	}
  
  plik>>tmps;
  if(tmps!="obszarow]" && tmps!="OBSZAROW]")
  { cerr<<"Brak znacznika [liczba obszarow] w msh"<<endl;  exit(1);	}
  else
   { 	plik>>tmpi;  	liczba_obszarow=tmpi; 		
 if (v==1) cout << "[liczba obszarow]\n"<< liczba_obszarow<<endl;
  }

  plik>>tmps;
  if(tmps!="[wspolrzedne" && tmps!="[WSPOLRZEDNE")
  { 	cerr<<"Brak znacznika [wspolrzedne wezlow] w msh"<<endl; exit(1);}

  plik>>tmps;
  if(tmps!="wezlow]" && tmps!="WEZLOW]")
  {	cerr<<"Brak znacznika [wspolrzedne wezlow]  w msh"<<endl; exit(1);}
  else
    {
    if (v==1) cout << "[wspolrzedne wezlow]"<<endl;
	double x,y; int ind;
	wspolrzedne=new wezel_[liczba_wezlow];
	for(int i=0;i<liczba_wezlow;i++)
	{
	  plik>>ind;	ind-=1;
	  plik>>x;	plik>>y;
	  wspolrzedne[ind].x=x;		wspolrzedne[ind].y=y;
 	if (v==1) cout << i+1 << "\t" << wspolrzedne[ind].x << "\t"
                       << wspolrzedne[ind].y<<endl;
	}
    }

  plik>>tmps;
  if(tmps!="[koneksje]" && tmps!="[KONEKSJE]")
  { cerr<<"Brak znacznika [koneksje] w msh"<<endl; exit(1);	}
  else
   {
     if (v==1) cout << "[koneksje]"<<endl;
	int tempik; 	koneksje=new koneksja_[liczba_elementow];
	for(int i=0;i<liczba_elementow;i++)
	 {
	  plik>>tempik;
	  koneksje[i].nr=tempik;
	  plik>>tmps;
	  if(tmps!="2D31")cerr<<"Blad koneksji -> element:"<<koneksje[i].nr<<endl;
	  koneksje[i].typ_el=tmps;
	  int tab[3];
	  plik>>tab[0];
	  plik>>tab[1];
	  plik>>tab[2];
	  plik>>tempik;
	  koneksje[i].obszar=tempik;
	  koneksje[i].wez_tw1=tab[0];
	  koneksje[i].wez_tw2=tab[1];
	  koneksje[i].wez_tw3=tab[2];
 if (v==1) cout << i+1 << "\t2D31\t"<<tab[0]<<"\t"<<tab[1]<<"\t"<<tab[2]<<"\t1" <<endl;
	}				
   }
		
  plik>>tmps;
  
  if(tmps!="[brzegi]" && tmps!="[BRZEGI]")
  { cerr<<"Brak znacznika [brzegi] w msh"<<endl;	exit(1);	}
  else
 {
 if (v==1) cout << "[brzegi]"<<endl;  
  brzegi=new brzeg_[liczba_brzegow];    int tempik;
	for(int i=0;i<liczba_brzegow;i++)
	   {
		plik>>tempik;
		brzegi[i].nr=tempik;
		plik>>tempik;
		brzegi[i].il_wez=tempik;
		plik>>tempik;
		brzegi[i].ob=tempik;
		brzegi[i].wt=new int [brzegi[i].il_wez];
 if (v==1) cout << brzegi[i].nr << " " << brzegi[i].il_wez << " " << brzegi[i].ob <<endl;

		for(int j=0;j<brzegi[i].il_wez;j++)
		{    plik>>tempik; 	brzegi[i].wt[j]=tempik; 
 		if (v==1)  cout << brzegi[i].wt[j] << " ";
                }
                if (v==1) cout << endl;
  
           }
	  
  }	  
	

plik>>tmps;
if(tmps!="[obszary]" && tmps!="[OBSZARY]")
{ cerr<<"Brak znacznika [OBSZARY] w msh)"<<endl;exit(1); }
else
{
  if (v==1) cout << "[obszary]" << endl;
  int tempik;
  obszary=new obszar_ [liczba_obszarow];
   for(int i=0;i<liczba_obszarow;i++)
	{
	  plik >> tempik;  obszary[i].nr=tempik;
	  plik >> tempik;  obszary[i].il_el=tempik;
          if (v==1) cout << obszary[i].nr << " " << obszary[i].il_el << endl;
	  //plik>>tmps;
	//  getline(plik,obszary[i].opis);
	  obszary[i].el_tw=new int [obszary[i].il_el];
	for (int j=0; j<obszary[i].il_el; j++)
           {   
              plik >> tempik;
              obszary[i].el_tw[j]=tempik;
              if (v==1) cout<<obszary[i].el_tw[j] << " ";
           }
             if (v==1) cout << endl;
         }
   }

   plik.close();
}
//============== Wczytywanie BC ===============
void Siatka::wczytaj_bc (string nazwa_pliku)
{
 nazwa_pliku+=".bc";
 ifstream plik; 
 plik.open (nazwa_pliku.c_str());
 
 string temp,temp2;
 
 plik >> temp; plik >> temp2;
 
 if (v==1)cout << "\nWczytywanie pliku bc: \n========================" << endl;

 if (temp+" "+temp2!="<WARUNKI BRZEGOWE>" && temp+" "+temp2!="<warunki brzegowe>") 
 { cerr << "Blad: nieprawidlowy format pliku bc (brak <WARUNKI BRZEGOWE>)"<<endl; exit(1); }
 else  {if (v==1)cout << "<WARUNKI BRZEGOWE> " <<endl; }

 while (temp!="ThElPlStressSimulator_1.0") 
 { 
   plik >> temp; 
   if (temp=="ThElPlStressSimulator_1.0" && v==1) cout << "\nThElPlStressSimulator_1.0\n"<<endl;
 }
 
 int tmp;
 plik >> tmp;
 if (tmp!=liczba_wymiarow) 
{ cerr << "Blad: Nieprawidlowa liczba wymiarow w pliku bc" << endl; exit(1); }
 else {if (v==1)cout << tmp << " "; }

 plik >> tmp;
 if (tmp!=liczba_wezlow)
 { cerr << "Blad: nieprawidlowa liczba wezlow w pliku bc" << endl; exit(1); }
 else { if (v==1)  cout << tmp<< " "; }
 
plik >> tmp;
 if (tmp!=liczba_elementow) 
{ cerr << "Blad: nieprawidlowa liczba elementow w pliku bc" << endl; exit(1); }
 else  {if (v==1)cout << tmp<< " ";}

 plik >> tmp;
 if (tmp!=liczba_obszarow) 
{ cerr << "Blad: nieprawidlowa liczba obszarow w pliku bc" << endl; exit(1); }
  else {if (v==1)cout << tmp << " ";}

 plik >> tmp;
 if (tmp!=liczba_brzegow) 
{ cerr << "Blad: nieprawidlowa liczba brzegow w pliku bc" << endl; exit(1); }
 else {if (v==1) cout <<tmp  << endl;}

 plik >> tmp; if (v==1) cout << tmp << endl;
 if (tmp!=liczba_brzegow)
  { cerr << "Blad: nieprawidlowa liczba brzegow w pliku bc" << endl; exit(1); }

 bc = new bc_ [liczba_brzegow];
 int idwar ,l_par;
 for (int i=0; i<liczba_brzegow; i++)
 {
   	plik >> bc[i].nr_brzegu; if (v==1) {cout << bc[i].nr_brzegu << " ";}
	plik >> bc[i].nr_obszaru; if (v==1) cout << bc[i].nr_obszaru << " ";
	plik >> bc[i].il_wezlow;  if (v==1)cout << bc[i].il_wezlow << " ";
	plik >> bc[i].il_warunkow; if (v==1)cout << bc[i].il_warunkow << endl;
 //---
for (int j=0; j<bc[i].il_warunkow; j++)
{
        plik >> tmp;
	if (tmp!=bc[i].il_wezlow)
	 {cerr << "Blad: nieprawidlowa liczba wezlow w pliku bc" << endl; exit(1);}
	 else if (v==1) cout << tmp << " ";
	plik >> idwar;

	if (v==1)cout << idwar << " "; 
        plik >> tmp;
	if (tmp!=0) { cerr << "Blad: nieprawidlowy parametr \"jak zdefiniowany\""<<endl;
	                exit(1); } else if (v==1)cout << tmp << " ";
			
	plik >> l_par; if (v==1)cout << l_par << endl;
        
	
	
if(idwar==1008)
{			
	   for (int k=0; k<l_par; k++)
	   {
	     plik >> idwar;
	     if (v==1)cout << idwar << " ";
	     int t;
	     plik >> t;  if (v==1) cout << t << " ";  plik >> t; if (v==1) cout << t << endl;
	     
	     if(idwar==1001) { 	double war;
				plik>>war;
				if (v==1)cout <<war << "\n" << endl;
				brzegi[bc[i].nr_brzegu].obc_x=true;
				brzegi[bc[i].nr_brzegu].wobc_x=war;
			     }
	     if(idwar==1002)
	            	     {
				double war;
				plik>>war;
				if (v==1)cout <<war << "\n" << endl;
				brzegi[bc[i].nr_brzegu].obc_y=true;
				brzegi[bc[i].nr_brzegu].wobc_y=war;
	                     }
	   //  else{cerr<<"Blad: dla warunku 1008 "<<endl; exit(1);}
	     
	   }		

}

else if(idwar==1002)
	{
	  brzegi[bc[i].nr_brzegu].utw_x=true;
	}

 

else if(idwar==1004)
	{
	  brzegi[bc[i].nr_brzegu].utw_y=true;
	}
else if(idwar==1006)
	{
	  brzegi[bc[i].nr_brzegu].utw_y=true;
	  brzegi[bc[i].nr_brzegu].utw_x=true;
	}
 else
	{
	  cerr<<"nieobslugiwany warunek brzegowy"<<endl;
	  exit(1);
	}
	

 }
 //---		
 }

}
//============== Wczytywanie M ================
void Siatka::wczytaj_m(string nazwa_pliku)
{
    nazwa_pliku+=".m";
 

  ifstream plik;
  plik.open (nazwa_pliku.c_str());
  string temp;
  if (v==1)cout << "\nWczytywanie pliku m \n==========================="<<endl;
  plik >> temp;
  
  if (temp!="<WLASNOSCI" && temp!="<wlasnosci")
  {
    cerr << "Blad: Nieprawidlowa forma znacznika <WLASNOSCI MATERIALOWE>" <<endl;
    exit (1);
  } else {if (v==1)cout << temp << " ";}
  
  plik >> temp;
  
  if (temp!="MATERIALOWE>" && temp!="materialowe>")
  {
       cerr << "Blad: Nieprawidlowa forma znacznika <WLASNOSCI MATERIALOWE>" <<endl;
    exit (1);
  } else {if (v==1) cout << temp << endl;}
 
  while (temp!="ThElPlStressSimulator_1.0") 
 { 
   plik >> temp; 
   if (temp=="ThElPlStressSimulator_1.0" && v==1) cout << "\nThElPlStressSimulator_1.0\n"<<endl;
 }
 
 int tmp;
 
 plik >> tmp;
 
 if (tmp!=liczba_wymiarow)
 {
   cerr << "Blad: nieprawidlowa liczba wymiarow w pliku m" << endl; exit(1);

 } else {if (v==1)cout << tmp << " ";}
 plik >> tmp;
  if (tmp!=liczba_wezlow)
 {
   cerr << "Blad: nieprawidlowa liczba wezlow w pliku m" << endl; exit(1);
 } else {if (v==1)cout << tmp << " ";}
 plik >> tmp;
  if (tmp!=liczba_elementow)
 {
   cerr << "Blad: nieprawidlowa liczba elementow w pliku m" << endl; exit(1);
 } else {if (v==1)cout << tmp << " ";}
 plik >> tmp;
  if (tmp!=liczba_obszarow)
 {
   cerr << "Blad: nieprawidlowa liczba obszarow w pliku m" << endl; exit(1);
 } else {if (v==1)cout << tmp << " ";}
 plik >> tmp;
  if (tmp!=liczba_brzegow)
 {
   cerr << "Blad: nieprawidlowa liczba brzegow w pliku m" << endl; exit(1);
 } else {if (v==1)cout << tmp << endl;}
  
  plik >> tmp;
	if(tmp!=701)
	{
		cerr<<"brak 701 w plik m"<< endl; exit(1);
	} else {if (v==1) cout << tmp << endl; }
  plik >>tmp;
	if(tmp!=5)
	{
		cerr<<"brak 5 w pliku m" << endl; exit(1);
	} else {if (v==1)cout << tmp << endl;}
  
  
  
  cout.precision(7);
  cout.flags(ios::scientific);
  
  for(int i=0;i<5;i++)
  { plik >> tmp;
    switch(tmp)
    {
     case 801: 
     { if (v==1)cout << "801 "; plik>>tmp; if (v==1){cout <<tmp<<endl;} plik >> w801; if (v==1)cout << w801<<endl; break;}
     case 802: 
     { if (v==1)cout << "802 "; plik>>tmp; if (v==1){cout <<tmp<<endl;} plik >> w802; if (v==1)cout << w802<<endl; break; }
     case 803: 
     { if (v==1)cout << "803 "; plik>>tmp; if (v==1){cout <<tmp<<endl;} plik >> w803; if (v==1)cout << w803<<endl; break; }
     case 804: 
     { if (v==1)cout << "804 "; plik>>tmp; if (v==1){cout <<tmp<<endl;} plik >> w804; if (v==1)cout << w804<<endl; break; }
     case 805: 
     { if (v==1)cout << "805 "; plik>>tmp; if (v==1){cout <<tmp<<endl;} plik >> w805; if (v==1)cout << w805<<endl; break; }
    }
    
  }

  
  plik.close();
}
//---------------------------------------------------------------------------------------------------
void Siatka::wspol_sprez(double e,double v)
{
  


  if(psn==1 && pso==0)
  {
	f1=(e)/(1-(v*v)); 
	f2=(e*v)/(1-(v*v));
	f3=(e)/(2*(1+v)); 
  }
  if(pso==1 && psn==0)
  {
	f1=(e*(1-v))/((1+v)*(1-(2*v))); 
	f2=(e*v)/((1+v)*(1-(2*v))); 
	f3=(e)/(2*(1+v));
  }
  if (v1==1) cout << "\n\nf1 = " << f1 << " f2 = " << f2 << " f3 = " << f3 << endl;
}
//==========================================================
double * Siatka:: wspolczynniki (int nr_elementu)
{
  double A; 
  double *tab=new double [7];
  double x1,y1, x2,y2, x3,y3;
  int tmpnw;
  
  tmpnw = koneksje[nr_elementu].wez_tw1;
  x1=wspolrzedne[tmpnw-1].x;
  y1=wspolrzedne[tmpnw-1].y;
  
  tmpnw = koneksje[nr_elementu].wez_tw2;
  x2=wspolrzedne[tmpnw-1].x;
  y2=wspolrzedne[tmpnw-1].y;
  
  tmpnw = koneksje[nr_elementu].wez_tw3;
  x3=wspolrzedne[tmpnw-1].x;
  y3=wspolrzedne[tmpnw-1].y;
  
  A=0.5*(x1*(y2-y3)+(x2*(y3-y1))+(x3*(y1-y2)));
  
  tab[0]=(y2-y3)/(2*A);//C21
  tab[1]=(y3-y1)/(2*A);//C22
  tab[2]=(y1-y2)/(2*A);//C23
  tab[3]=(x3-x2)/(2*A);//C31
  tab[4]=(x1-x3)/(2*A);//C32	
  tab[5]=(x2-x1)/(2*A);//C33
  tab[6]=A;//pole elementu

  
return tab;	
}
//===========================================================

void Siatka::Macierz()
{
  int rozmiar=liczba_wezlow*2;
  //!!!!!!!!!!!!!!!!  ++2
  mac =new double* [rozmiar+2];
  for (int i=0; i<rozmiar+1; i++) mac[i] = new double [rozmiar+2];
  
  //zeruje nowa macierz
  
  for (int i=0; i<rozmiar; i++)
  for (int j=0; j<rozmiar; j++) mac[i][j]=0;
  
  for(int i=0;i <liczba_elementow; i++)
  {

    int nr_el=koneksje[i].nr-1;
    double *wsp=wspolczynniki(nr_el);
    int tab_wez[6];
    tab_wez[0]=(koneksje[nr_el].wez_tw1*2)-2;//i11
    tab_wez[1]=(koneksje[nr_el].wez_tw1*2)-1;//i12
    tab_wez[2]=(koneksje[nr_el].wez_tw2*2)-2;//i21
    tab_wez[3]=(koneksje[nr_el].wez_tw2*2)-1;//i22	
    tab_wez[4]=(koneksje[nr_el].wez_tw3*2)-2;//i31	
    tab_wez[5]=(koneksje[nr_el].wez_tw3*2)-1;//i12
    
    
  
int a=0;

for(int x=0;x<3;x++)
{
int b=0;
for(int z=0;z<3;z++)
{				
mac[tab_wez[a]][tab_wez[b]]+=wsp[6]*((f1*wsp[x]*wsp[z])+(f3*wsp[3+x]*wsp[3+z]));	
mac[tab_wez[a]][tab_wez[b+1]]+= wsp[6]*((f2*wsp[x]*wsp[3+z])+(f3*wsp[3+x]*wsp[z]));			
mac[tab_wez[a+1]][tab_wez[b]]+=wsp[6]*((f2*wsp[3+x]*wsp[z])+(f3*wsp[x]*wsp[3+z]));
mac[tab_wez[a+1]][tab_wez[b+1]]+= wsp[6]*((f1*wsp[3+x]*wsp[3+z])+(f3*wsp[x]*wsp[z]));

b+=2;
}
a+=2;
}		
delete [] wsp;

    
  }

if (m==1)
{
  ofstream plik;
  plik.precision(7);
  plik.flags(ios::scientific);
  plik.open (nazwam.c_str());
  for (int i=0; i<rozmiar; i++)
    {
        for (int j=0; j<rozmiar; j++)
         {
           plik << mac[i][j]<<"\t";
         } 
      plik << endl;
    }
  cout << "\n --> macierz sztywnosci zostala zapisana do pliku: " << nazwam << endl;
  plik.close();
}

}
//========================================================
void Siatka::Obciazenia ()
{
  //zeruje wektor;
 for (int i=0; i<liczba_wezlow*2; i++) mac[i][(liczba_wezlow*2)]=0;
 //---------------------------------------------
   
   for (int i=0; i<liczba_brzegow; i++)
   {
      if(brzegi[i].obc_x==true) 
      {
        for(int j=0;j<(brzegi[i].il_wez-1);j++)
	  {
            double h,x1,x2,y1,y2; 
            
             x1=wspolrzedne[brzegi[i].wt[j]-1].x;
             x2=wspolrzedne[brzegi[i].wt[j+1]-1].x;
   
             y1=wspolrzedne[brzegi[i].wt[j]-1].y;
             y2=wspolrzedne[brzegi[i].wt[j+1]-1].y;

             h=sqrt(((x1-x2)*(x1-x2))+((y1-y2)*(y1-y2)));

 	      mac[((brzegi[i].wt[j])*2)-2][(liczba_wezlow*2)]+=(h/2)*brzegi[i].wobc_x;
              mac[((brzegi[i].wt[j+1])*2)-2][(liczba_wezlow*2)]+=(h/2)*brzegi[i].wobc_x;
          }
      }

      if(brzegi[i].obc_y==true) 
      {
        for(int j=0;j<(brzegi[i].il_wez-1);j++)
	  {
            double h,x1,x2,y1,y2; //,q2,q1;

             x1=wspolrzedne[brzegi[i].wt[j]-1].x;
             x2=wspolrzedne[brzegi[i].wt[j+1]-1].x;
   
             y1=wspolrzedne[brzegi[i].wt[j]-1].y;
             y2=wspolrzedne[brzegi[i].wt[j+1]-1].y;

             h=sqrt(((x1-x2)*(x1-x2))+((y1-y2)*(y1-y2)));

		mac[((brzegi[i].wt[j])*2)-1][(liczba_wezlow*2)]+=(h/2)*brzegi[i].wobc_y;
		mac[((brzegi[i].wt[j+1])*2)-1][(liczba_wezlow*2)]+=(h/2)*brzegi[i].wobc_y;

          }
      }
   }
if (v1==1)
{
  cout << "Wektor obciazen:"<<endl;
  for (int i=0; i<liczba_wezlow*2; i++) cout << i <<"  "<< mac[i][liczba_wezlow*2] << endl;

}


}
//==================================================
void Siatka::Utwierdzenia()
{

 for(int i=0;i<liczba_brzegow;i++)
 {
  if(brzegi[i].utw_x)
  {
  for(int j=0;j<brzegi[i].il_wez;j++)
  {
    int ind=(brzegi[i].wt[j]*2)-2;
    mac[ind][liczba_wezlow*2]=0;
    mac[ind][ind]=1;
    for(int k=(ind+1);k<(liczba_wezlow*2);k++) mac[ind][k]=0;
    for(int k=(ind-1);k>=0;k--) mac[ind][k]=0;
    for(int k=(ind+1);k<(liczba_wezlow*2);k++) mac[k][ind]=0;
    for(int k=(ind-1);k>=0;k--) mac [k][ind]=0;
  }
  }
  if(brzegi[i].utw_y)
  {
  for(int j=0;j<brzegi[i].il_wez;j++)
  {
    int ind=(brzegi[i].wt[j]*2)-1;
    mac[ind][liczba_wezlow*2]=0;
    mac[ind][ind]=1;
    for(int k=(ind+1);k<(liczba_wezlow*2);k++) mac[ind][k]=0;
    for(int k=(ind-1);k>=0;k--) mac[ind][k]=0;
    for(int k=(ind+1);k<(liczba_wezlow*2);k++) mac[k][ind]=0;
    for(int k=(ind-1);k>=0;k--) mac [k][ind]=0;
  }
 
  }
 }

if (v1==1){
 cout << "Wektor obciazen po utwierdzeniu:"<<endl;
  for (int i=0; i<(liczba_wezlow*2); i++)
    cout << mac[i][liczba_wezlow*2] <<endl;
}

}
//===========================================================

double * Siatka::Gauss ()
{

  double suma;
  int i, j, k, rozmiar=liczba_wezlow*2;
  double *tab;
  tab = new double [rozmiar];

  // ELIMINACJA ZMIENNYCH - tworzenie macierzy trojkatnej gornej

  for (k=0; k<rozmiar-1; k++)
  {
     for (i=k+1; i<rozmiar; i++)
       for (j=rozmiar; j>=k; j--)
          mac[i][j] = mac[i][j] - mac[i][k]*mac[k][j]/mac[k][k];
  }

  // OBLICZANIE WEKTORA ROZWIAZAN

  for (i=rozmiar-1; i>=0; i--)
  {
    suma=0.0;
    for (j=rozmiar-1; j>=i+1; j--)
     suma = suma + mac[i][j]*mac[j][rozmiar];
     mac[i][rozmiar] = (mac[i][rozmiar]-suma)/mac[i][i];

  }  
if (v==1)
{
 cout << "Przemieszczenia: "<<endl;
 for (int i=0; i<rozmiar; i=i+2)
 {
 
  cout << mac[i][rozmiar] << "\t" << mac[i+1][rozmiar]<<endl;
 }
}
 for (int i=0; i<rozmiar; i++)
 {
 
  tab[i] =  mac[i][rozmiar];
 }


return tab;

/*

	double * temp=new double[liczba_wezlow*2];
	for(int i=0;i<(liczba_wezlow*2);i++)temp[i]=0;
	for(int i=1;i<(liczba_wezlow*2);i++)
	{
               for(int k=i;k<(liczba_wezlow*2);k++)
               {
               	for(int j=i;j<((liczba_wezlow*2)+1);j++) mac[k][j]=mac[k][j]-(mac[i-1][j]*(mac[k][i-1]/mac[i-1][i-1]));
               	mac[k][i-1]=0;
               }
		for(int i=((liczba_wezlow*2)-1);i>=0;i--)
		{
		       double tmp=0;
		       for(int k=i+1;k<(liczba_wezlow*2);k++) tmp+=mac[i][k]*temp[k];
               		temp[i]=(mac[i][liczba_wezlow*2]-tmp)/mac[i][i];
		}
	}
return temp;*/


}
//====================================================================
odksztalcenie_ * Siatka::Odksztalcenia (double v)
{
odksztalcenie_ *tab;
   tab = new odksztalcenie_ [liczba_elementow];

   if(psn==1) 
    {
       for (int i=0; i<liczba_elementow; i++)
        {
            int nr_el=koneksje[i].nr-1;

   
         	 
 		 int wt1,wt2,wt3;
                 double u1,u2,u3,v1,v2,v3;

                 wt1=koneksje[i].wez_tw1;
                 wt2=koneksje[i].wez_tw2;
                 wt3=koneksje[i].wez_tw3;
                 int k = (liczba_wezlow*2);

		u1 = mac[(wt1*2)-2][k];
     		u2 = mac[(wt2*2)-2][k];
 		u3 = mac[(wt3*2)-2][k];
 	
		v1 = mac[(wt1*2)-1][k];
 		v2 = mac[(wt2*2)-1][k];
 		v3 = mac[(wt3*2)-1][k];
               
                double *tabWSP=wspolczynniki(nr_el);

  tab[i].Ex=(tabWSP[0]*u1)+(tabWSP[1]*u2)+(tabWSP[2]*u3);
  tab[i].Ey=(tabWSP[3]*v1)+(tabWSP[4]*v2)+(tabWSP[5]*v3);
  tab[i].Qxy=(tabWSP[3]*u1)+(tabWSP[4]*u2)+(tabWSP[5]*u3)+(tabWSP[0]*v1)+(tabWSP[1]*v2)+(tabWSP[2]*v3);
  tab[i].Ez=((-1*v)/(1-v))*(tab[i].Ex+tab[i].Ey);
                delete []tabWSP;
	
        }

    }

   else
    {
       for (int i=0; i<liczba_elementow; i++)
          {
             int nr_el=koneksje[i].nr-1;

 		 int wt1,wt2,wt3;
                 double u1,u2,u3,v1,v2,v3;

                 wt1=koneksje[i].wez_tw1;
                 wt2=koneksje[i].wez_tw2;
                 wt3=koneksje[i].wez_tw3;
                 int k = (liczba_wezlow*2);
 		
		u1 = mac[(wt1*2)-2][k];
     		u2 = mac[(wt2*2)-2][k];
 		u3 = mac[(wt3*2)-2][k];
 	
		v1 = mac[(wt1*2)-1][k];
 		v2 = mac[(wt2*2)-1][k];
 		v3 = mac[(wt3*2)-1][k];
               
                double *tabWSP=wspolczynniki(nr_el);

  tab[i].Ex=(tabWSP[0]*u1)+(tabWSP[1]*u2)+(tabWSP[2]*u3);
  tab[i].Ey=(tabWSP[3]*v1)+(tabWSP[4]*v2)+(tabWSP[5]*v3);
  tab[i].Qxy=(tabWSP[3]*u1)+(tabWSP[4]*u2)+(tabWSP[5]*u3)+(tabWSP[0]*v1)+(tabWSP[1]*v2)+(tabWSP[2]*v3);
  tab[i].Ez=0;

   		delete []tabWSP;
          }
    }
    
 return tab;
}
//======================================
naprezenie_ *Siatka::Naprezenia (double v, odksztalcenie_ *odk)
{
  naprezenie_ *tab;
  tab=new naprezenie_ [liczba_elementow];

   if(pso) 
   { 
     for (int i=0; i<liczba_elementow; i++)
        {

                tab[i].Sx=(f1*odk[i].Ex)+(f2*odk[i].Ey);
                tab[i].Sy=(f2*odk[i].Ex)+(f1*odk[i].Ey);
                tab[i].Yxy=f3*odk[i].Qxy;
                tab[i].Sz=v*(tab[i].Sx+tab[i].Sy);
        }
   }

   else
   {
      for (int i=0; i<liczba_elementow; i++)
        {

                tab[i].Sx=(f1*odk[i].Ex)+(f2*odk[i].Ey);
                tab[i].Sy=(f2*odk[i].Ex)+(f1*odk[i].Ey);
                tab[i].Yxy=f3*odk[i].Qxy;
                tab[i].Sz=v*0;
        }     
   }


  return tab;
}
//---------------------------------------------------------------------------
void Siatka::ZapiszOdksztalcenia (odksztalcenie_ * o)
{
  ofstream plik;
  plik.precision(7);
  plik.flags(ios::scientific);
  plik.open (nazwao.c_str());
plik << "<WYNIKI>\n\n[opis]\n\n[liczba wymiarow]\n"<<liczba_wymiarow<<"\n\n[liczba wezlow]\n"
     << liczba_wezlow << "\n\n[liczba elementow]\n"<<liczba_elementow<<"\n\n[wartosci w elementach]\n1\n\n"
     << "[liczba skladowych]\n4\n\n[nazwa wielkosci]\nodksztalcenia ex, ey, ez ,qxy\n\n[wartosci]"<<endl;

for (int i=0; i<liczba_elementow; i++)
{
  plik << i+1 <<"\t\t"<< o[i].Ex << "\t\t" << o[i].Ey << "\t\t" << o[i].Ez << "\t\t" << o[i].Qxy  << endl;
}


  plik.close();
}
void Siatka::ZapiszNaprezenia (naprezenie_ * n)
{
  ofstream plik;
  plik.precision(7);
  plik.flags(ios::scientific);
  plik.open (nazwan.c_str());
plik << "<WYNIKI>\n\n[opis]\n\n[liczba wymiarow]\n"<<liczba_wymiarow<<"\n\n[liczba wezlow]\n"
     << liczba_wezlow << "\n\n[liczba elementow]\n"<<liczba_elementow<<"\n\n[wartosci w elementach]\n1\n\n"
     << "[liczba skladowych]\n4\n\n[nazwa wielkosci]\nodksztalcenia sx, sy, sz, txy [Pa]\n\n[wartosci]"<<endl;

for (int i=0; i<liczba_elementow; i++)
{
  plik << i+1 <<"\t\t"<< n[i].Sx << "\t\t" << n[i].Sy << "\t\t" << n[i].Sz << "\t\t" << n[i].Yxy  << endl;
}


  plik.close();
}
void Siatka::ZapiszPrzemieszczenia (double *tab)
{
  ofstream prz;
  prz.precision(7);
  prz.flags(ios::scientific);
  prz.open (nazwap.c_str());

  prz << "<WYNIKI>\n\n[opis]\n\n[liczba wymiarow]\n2\n\n[liczba wezlow]\n" 
      << liczba_wezlow << "\n\n[liczba elementow]\n"<<liczba_elementow <<"\n\n[wartosci w elementach]\n0"
      << "\n\nliczba skladowych]\n2\n\n[nazwa wielkosci]\nu,v [m]\n\n[wartosci]" << endl;

  int licznik=0;
  for (int i=0; i<(liczba_wezlow*2); i=i+2)
   {
     prz << licznik+1 << "\t\t" << tab[i] << "\t\t" << tab[i+1]<<endl;
     licznik ++;
   }

  prz.close();
}
//---------------------------------------------------------------------------------------------------

//----------------------------------------------------------------------------------------------------
int main (int argc, char ** argv)
{
 
 Siatka s;
for (int i=0; i<argc; i++)
{
  string temp=argv[i];
  if (temp=="-v") { v=1; v1=1;}
  else if (temp=="-m") { m=1; nazwam=argv[i+1]; }
  else if (temp=="-psn") { s.setPSN(); }
  else if (temp=="-pso") { s.setPSO(); }
  else if (temp=="-all") { nazwa=argv[i+1]; }
  else if (temp=="-o") { 
                    string temp2=argv[i+1];
                    nazwan=temp2+"_n.rlt";
                    nazwap=temp2+"_p.rlt";
                    nazwao=temp2+"_o.rlt";
                   }
  else if (temp=="-on") {nazwan=argv[i+1];}
  else if (temp=="-op") {nazwap=argv[i+1];}
  else if (temp=="-oo") {nazwao=argv[i+1];}
  else if (temp=="-help")
    {
  
      cout << "Programowanie Aplikacji Inzynierskich v0.0.9 Beta\n"
           << "                                    by Maciej Kunowski\n\n"
           << "                    POMOC                        \n\n"
           << " parametry wywolania programu:\n\n"
           << " -v \t\t\t Wyswietla czynnosci, ktore wykonuje aplikacja\n"
           << " -m nazwa_pliku \t Zapisuje macierz sztywnosci do \"nazwa_pliku\"\n"
           << " -psn \t\t\t Ustawia plaski stan naprezen\n"
           << " -pso \t\t\t Ustawia plaski stan odksztalcen\n"
           << " -all nazwa_pliku \t Nazwa dla wszystkich plikow wejsciowych (bez rozszerzenia)\n"
           << " -o nazwa_pliku \t Nazwa dla wszystkich plikow wyjsciowych (bez rozszerzenia)\n"
           << " -on nazwa_pliku \t Nazwa dla pliku naprezen (z rozszerzeniem)\n"
           << " -op nazwa_pliku \t Nazwa dla pliku przemieszczen (z rozszerzeniem)\n"
           << " -oo nazwa_pliku \t Nazwa dla pliku odszksztalcen (z rozszerzeniem)\n"<<endl;
         exit(1);
    }
 else if (temp[0]!='-') {
                           int d = temp.length();
                           if (temp[d-1]=='h' && temp[d-4]=='.') { 
							            temp.erase (d-4,d-1);
                                                                    nazwa_msh=temp;
                                                                 }
                           if (temp[d-1]=='c' && temp[d-3]=='.') {  
								    temp.erase (d-3,d-1);
									nazwa_bc=temp;
                                                                 }
			   if (temp[d-1]=='m' && temp[d-2]=='.') { 
                 						    temp.erase (d-2,d-1);
									nazwa_m=temp;
                                                                 }                                 
                        }
}


if (nazwa_msh!="pusty") s.wczytaj_msh (nazwa_msh.c_str());
else s.wczytaj_msh (nazwa.c_str());

if (nazwa_bc!="pusty") s.wczytaj_bc (nazwa_bc.c_str());
else s.wczytaj_bc (nazwa.c_str());

if (nazwa_m!="pusty") s.wczytaj_m (nazwa_m.c_str());
else s.wczytaj_m (nazwa.c_str());

/*
if (nazwa!="pusty")
{
  s.wczytaj_msh(nazwa.c_str());
  s.wczytaj_bc (nazwa.c_str());
  s.wczytaj_m  (nazwa.c_str()); 
}
else
{
  cout << "Podaj nazwe pliku *.msh (bez rozszerzenia): ";
  cin >> nazwa_msh;
  cout << "Podaj nazwe pliku *.bc (bez rozszerzenia): ";
  cin >> nazwa_bc;
  cout << "Podaj nazwe pliku *.m (bez rozszerzenia): ";
  cin >> nazwa_m;
  s.wczytaj_msh(nazwa_msh.c_str());
  s.wczytaj_bc(nazwa_bc.c_str());
  s.wczytaj_m(nazwa_m.c_str());
}
*/
 s.wspol_sprez (s.w_801_(), s.w_804_());
 s.Macierz();
 s.Obciazenia();
 s.Utwierdzenia();
 //----------
 double *wek = s.Gauss();
 odksztalcenie_ *odk = s.Odksztalcenia (s.w_804_());
 naprezenie_ *nap = s.Naprezenia (s.w_804_(), odk);
//-----------
 s.ZapiszOdksztalcenia(odk);
 s.ZapiszNaprezenia(nap);
 s.ZapiszPrzemieszczenia(wek);
 return 0;
}
